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Abstract 

This paper presents an approach to modeling the dynamics of flexible multibody systems such 
as flexible spacecraft and limber space robotic systems. A large number of degrees of freedom 
and complex dynamic interactions are typical in these systems. This paper uses spatial operators 
to develop efficient recursive algorithms for the dynamics of these systems. This approach very 
efficiently manages complexity by means of a hierarchy of mathematical operations. 

1. Introduction 

A wide range of complex mechanical systems can be modeled as a set of hinge-connected flexible 
and rigid bodies. This paper presents an approach to modeling the dynamics of such systems that 
uses spatial operators. This approach very efficiently manages complexity by means of a hierarchy 
of mathematical operations. The highest level in this hierarchy consists of spatial operators which 
relate velocities, accelerations, and forces between distinct points in the system. At lower levels, each 
spatial operator is decomposed easily into detailed spatially recursive algorithms to do computation. 
The recursive algorithms are cast within the highly developed framework of filtering and smoothing 
theory. Algorithms which are quite popular in state estimation theory for discrete-time systems 
can now be applied to conduct spatially recursive operations essential in multibody dynamics. The 
main focus is on serial chains, but extensions to general topologies are also described. Comparison 
of computational costs illustrates the efficiency of the recursive algorithms. 

This paper uses spatial operators [1,2] to develop efficient recursive algorithms for flexible multibody 
systems for such applications as flexible spacecraft and limber space robotic systems. A large 
number of degrees of freedom and complex dynamic interactions are typical in these systems. The 
main contributions of the paper are: (1) high-level architectural understanding of the mass matrix 
and its inverse, (2) high-level expressions which can be easily implemented with spatial Kalman 
filtering and smoothing, (3) efficient inverse and forward dynamics recursive algorithms, and (4) 
analysis of computational cost of the new algorithms. This adds to the rapidly developing body of 
research in the recursive dynamics of flexible multibody systems [3-5]. 


2. Equations of Motion 

Equations of motion are developed for a serial system formed by N articulated flexible bodies. 
Recursive relationships between the modal velocities, accelerations and forces are developed. Spatial 
operators express these relationships compactly to obtain what is referred to here as the Newton- 
Euler mass matrix factorization. 

Each flexible body has a lumped mass model formed by a set of nodal rigid bodies. Such models 
are typically developed using standard finite element analysis. The k th body has n 3 (k) nodes. The 
j th node on the k th body is called the jj! 1 node. There is a body reference frame Tk for the k th body. 
Deformation of the nodes on the body is described with respect to this body reference frame, while 
the rigid-body motion of the k th body is characterized by the motion of frame Tk- 

The 6-dimensional spatial deformation (slope plus translational) of node jk (with respect to frame 
jF fc ) is u(jk ) € ft 6 . The overall deformation field for the k th body is the vector u(k) = col|u(jfc)| € 
The vector from frame Tk to the reference frame on node jk is l{k,jk ) G ft 3 - 
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The spatial inertia of the j th node is 


M s (j k ) = 


J(jk) m (jk)p(jk) 

- ™(jk)P(jk ) rn{j k )I 


€ ft 6x6 


( 1 ) 


where J(jk ), p(jk) and m(jk) are the inertia tensor about the node reference frame, the vector 
from the node reference frame to its center of mass, and the mass, respectively, for the j th node on 
the k th body. The structural mass matrix for the k th body M s (k ) is the block diagnal matrix 

M s {k) = diag {M s {jk)} € » 6n *«* 6n ‘W (2) 

The structural stiffness matrix is denoted K s (k) € & 6n *( fc ) x6n »(*>. For a 3-vector x, there is a 
corresponding cross product matrix x. Both M s (k) and K 3 ( fc ) are typically generated using finite 
element analysis. 

As in Figure 1, the bodies in the serial chain are numbered in increasing order from tip to base. The 



Towards base -« ► Towards tip 

Figure 1: Illustration of links and hinges in a flexible serial multi- 
body system 

terms inboard (outboard) denote the direction along the serial chain towards (away from) the base 
body. The k th body is attached on the inboard side to the ( k + l) th body by the k th hinge, and on 
the outboard side to the ( k - l) th body by the (k - l) tb hinge. On the k th body, the node to which 
the outboard hinge (the ( k - l) ih hinge) is attached is node t k , while the node to which the inboard 
hinge (the k th hinge) is attached is node d k . The k th hinge couples nodes d k and t*+i- Attached 
to each of these nodes are the k th hinge reference frames O k and O £ respectively. The number of 
degrees of freedom (dofs) for the k th hinge is n r (k). The vector of configuration variables for the 
k th hinge is 6(k) 6 ^Sl nr< ' k \ while the hinge’s vector of generalized speeds is (3(k ) e & nr ^. In general, 
when there axe nonholonomic hinge constraints, the dimensionality of /3 (k) may be less than that of 
0(k). For convenience, and without any loss in generality, it is assumed here that the dimensions of 
the vectors 6(k) and j3(k) are equal. In most situations, fi(k) is simply 6. However there are many 
cases where the use of quasi-coordinates simplifies the dynamical equations of motion, and there 
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may be a better choice for /3(fc). The relative spatial velocity A v (k) across the hinge is H*(k)/3(k), 
where H*(k ) is the joint map matrix for the k th hinge. 

A set of n m {k) assumed modes is chosen for the k th body. Let n^(fc) € 3? 6 be the modal spatial 
displacement vector at the f k h node for the r th mode. The modal spatial displacement influence 
vector IF(fc) e 3? 6xn ™W for the jjj 1 node and the modal matrix II(A;) € & 6n 4*0><"m( fc ) for the 
k th body are 

n j (k) = [ni(fc), • • , ni m(Jb) (*)] and n(*) = coi{ip'(fc)} 

The r th column of II(Jt) is II r (A:), which defines the mode shape for the r th assumed mode for the 
k th body. Let 77(A) € be the vector of modal deformation variables for the k th body. The 

spatial deformation of node j k and the spatial deformation field u(Ar) for the A;* body are 

u(jk) = IP {k)r](k) and u(k ) = n(k)r/(k) (3) 

Note that for cantilever modes 

II?(Jfc) = 0 and r=l---n m (A:) (4) 

The vector of generalized configuration variables i?( k ) and the vector of generalized speeds x{k) for 
the k th body are 

m = ( e xw = ( ) s « 

where Af(k) = n m (k ) + n T (k). The overall vectors of generalized configuration variables 0 and 
generalized speeds x f° r the serial multibody system are 

d = col{t?(Jb)} 6 3^ and x = col{x(A:)} e (6) 

where M = 52k =1 A/” (A;). The number of overall degrees of freedom for the multibody system is 
A f. The state of the multibody system is defined by the pair of vectors {t?,.x}- For a given system 
state {i?,x}> the equations of motion relate the generalized accelerations x and generalized forces 
T € fR*. The inverse dynamics problem is to compute the generalized forces T for a prescribed 
set of generalized accelerations x ■ Conversely, the forward dynamics problem is to compute the 
generalized accelerations x from the generalized forces T. 

2.1 Recursive Propagation of Velocities 

Let V(k) be the spatial velocity of the k th body reference frame At, i.e., 

where w(Ar) and v(k) are the angular and linear velocities respectively of T k - The spatial velocity 
V s (t k + 1) € 3£ 6 of node t k + 1 (on the inboard side of the k th hinge) is related to the spatial velocity 
V(k + 1) of the (k + 1 ) tfl body reference frame T k + 1, and the modal deformation variable rates 
V(k+1): 

V,(t k+ 1) = f(HM k+ i)nHl) + n'(Hl)ij(Hl) ( 7 ) 
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( 8 ) 


The spatial transformation operator 4>(x,y) £ S? 6x6 is 

«*■»>-(£ V) 

where l(x,y ) £ 3? 3 is the vector between the points x and y. Note the group property 

4>{x,y)4>{y,z) = 4>{x,z) 

for arbitrary points x, y and z. As in Eq. (7), and all through this paper, the index k will be used 
to refer to both the k th body as well as to the k th body reference frame T k with the specific usage 
coming from the context. For instance, V(k) and <t>(k,tk ) are the same as V(T k ) and <f>(F k ,tk) 
respectively. 

The spatial velocity V(0 % ) of frame O* (on the inboard side of the k ih hinge) is related to V s (/*+i) 
by 

V(Ot) = <f(t k+ i,O k )V s (t k+1 ) (9) 

The relative spatial velocity Ay(fc) across the k th hinge is H m (k)p(k), and so the spatial velocity 
V(O k ) of frame O k on the outboard side of the k th hinge is 

V(O k ) = V(Ot) + H*(k)f3(k ) (10) 

The spatial velocity V(k ) of the k th body reference frame is 

V(k ) = 4>*(O k ,k)V(Ok) - u(d k ) = 4 >*(o k ,k) \v(Ok) - ir'wft*)] ( 11 ) 

Eq. (7), Eq. (9), Eq. (10) and Eq. (11) together imply 

v(k) = ^(Jk + i^Mfc + ^ + ^^fc+x^jn^ + i^ + i) 

+<f>*(Ok,k ) [jy*(*)/3(fc) - n d (k)f,(k)} (i2) 

Thus, with W{k) = n m (k ) + 6 and Eq. (12), the modal spatial velocity V m (k) £ for the 

k th body is 

v m(k)= ( ) = $*(k + i,k)v m (k + i) + n*(k)x(k)£sFM (13) 

where the interbody transformation operator $(., .) and the modal joint map matrix Ti(k) are 

*(* + !,*) S (14) 

e ^(k)xW(k) 


n(k) * f 1 -wm- 

n[k) - [ 0 Hjr(k) 


(15) 


where 

Note that 
where 


Hjr(k) = H(k)<t>(O k ,k)£U n 'W* fi , and Il£(A:) = <t>*(O k ,k)E d (k) £ » 6xAr « 

$(k + l,k) = A(k + l)B(k+l,k) (16) 


A(k) = ( [ f(k% ) G ^ (fc)X6 and B{k +1 ' k) = #**+!,*)] 6 ( 17 > 
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Also, the modal joint map matrix H(k) can be partitioned as 

«(*) = ( 


( 18 ) 


where 


«,(*) = [/, -[njwr] € and «,(*) = [«, (19) 

With AT = DfcLi the spatial operator £$ is defined as 


£* = 


/ 0 0 0 

$(2,1) 0 

0 $(3,2) ... 


V o 


0 


$ = [/ — £$] 1 = I + £* + • • • + £$ 1 — 


0 

0 ^ 


0 

0 


0 

0 

€ 

(N,N-1) 

o J 


le spatial operator $ as 

/ 

0 

... 0 \ 

$(2,1) 

I 

... 0 

$(AM) 

$(A,2) ... i) 


( 20 ) 


jjA f xJJ 


( 21 ) 


where 


$(i,j) i $(*,*- 1 ) ••• $(j + l,j) for t>i 


( 22 ) 


Also define the spatial operator H = diag|?f(fc)| G It follows that 

= W X 

where V m = coljV m (/c) j G 9^ • 

2.2 Modal Mass Matrix for a Single Body 

An expression for the modal mass matrix of the k th body is derived. Denote by V a (j k ) € £ 6 the 
spatial velocity of node j k . V s {k) = col{F s (j fc )} G % 6n ’ ik) is the vector of all nodal spatial velocities 

(23) 


for the k ih body. It follows (as in Eq. (7)) that 


V s (k) = B*(k)V( k) + u(k) = P(fc), B'(k))V m {k) 


where 


B(k) = [4>(k,l k ),<t>{k,2 k ),---,4>{k,n s {k))] G (24) 

Since M a (k) is the structural mass matrix of the k th body, the kinetic energy of the k th body is 


±V;(k)M 3 (k)V s (k) = \v^(k)M m (k)V m (k) 


where 


M m (k) = 


n m ) «■(*)] = ( m§ ( tj lew ) € RAWxW) <25) 
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Corresponding to the generalized speed vector x(^)» M m (k) is the modal mass matrix of the 
k th body. In the block partitioning in Eq. (25), the superscripts / and r denote the flexible and 
rigid blocks respectively. Thus M (k) represents the flex/flex coupling block, while M^(k) the 
flex/rigid coupling block, of M m {k). Note that M£(k) is precisely the rigid body spatial inertia of 
the k th body. Indeed, M m (k ) reduces to the rigid body spatial inertia when the body flexibility is 
ignored, i.e., no modes are used, since in this case n m (k) = 0 (and !!(&) is null). 

Since the vector l(k,jk ) from Tk to node jk depends on the deformation of the node, the operator 
B{k ) is also deformation dependent. From Eq. (25) it follows that while the block TiK ( k ) is 
deformation independent, both the blocks M^{k) and M£(k) are deformation dependent. The 
detailed expression for the modal mass matrix can be defined using modal integrals which are 
computed as a part of the finite-element structural analysis of the flexible bodies. These expressions 
for the modal integrals and the modal mass matrix of the k th body can be found in [6]. Often the 
deformation dependent parts of the modal mass matrix are ignored, and free-free eigen-modes are 
used for the assumed modes II(k). When this is the case, M^[(k) is zero and M^/(k) is block 
diagonal. 


2.3 Recursive Propagation of Accelerations 

Differentiation of the velocity recursion equation, Eq. (13), results in the following recursive ex- 
pression for the modal spatial acceleration a m (k ) 6 f or the k th body: 


a m (k) = V m (k) 


v( k ) 

a(k) 


= $*(k + 1 ,k)a m (k + 1) + 7-T(fc)x(k) + a m (*0 


where a(k ) = V(k), and the Coriolis and centrifugal acceleration term a m (k) £ is 


®m(k) — 


d$*(fc + l,k) 


V m (k + 1) + 


dH‘(k) 


X(k) 


(26) 


(27) 


d t '"‘ v " ' ' d t 

The detailed expressions for a m (k ) can be found in [6]. Defining a m = col|a OT (fc)j £ and 
a m = col|a m (fc)| £ Eq. (26) can be reexpressed using spatial operators in the form 

= {fli X "b ®m) (28) 

The vector of spatial accelerations of all the nodes for the k th body, a s (k) = col|o s (jA:)| 

6 is obtained by differentiating Eq. (23): 


where 


<*,(*) = V s (k) = [H(*), B*(k)]a m (k) + a(k) 

a(k ) = col{a(j fc )} = e 


(29) 

(30) 


2.4 Recursive Propagation of Forces 

The equations of motion for the k th body are now developed. Let f(k — 1) £ 3? 6 denote the effective 
spatial force of interaction, referred to frame Tk- 1 , between the k th and ( k — l) £/l bodies across the 
(k - l) £/l hinge. Recall that the ( k - l) £/l hinge is between node 4 on the k th body and node dk-i 
on the ( k — l) £/l body. With f s {jk) £ 3ft 6 denoting the spatial force at a node jk, the force balance 
equation for node 4 is 

/»(4) = ^(4,fc- l)/(k- l) + M s (4)a s (4) + 6(4) + /A-(4) (31) 
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For all nodes other than node t k on the k ih body, the force balance equation is 

fs(jk) = M s (j k )a s (j k ) + b(j k ) + fh'(jk) ( 32 ) 

In the above, f K {k) = K s (k)u(k ) € is the vector of spatial elastic strain forces for the nodes 

on the k th body, while b(j k ) € ft 6 is the spatial gyroscopic force for node j k and is given by 


v(jk)J(jk)u(j k ) 

m(jk MjkMjkWk) 


where v(jk) G ft 3 denotes the angular velocity of node jk- Define 

0 > 

4>(tk,k- 1) G & 6n *W x6 and b(k) = col{b(;*)} € ft 6ns(fc) (34) 

0 j 

Eq. (31) and Eq. (32) imply 

f s {k ) = C(k , k - 1 )f(k - 1) + M 3 (k)a a (k) + b(k ) + K s (k)u(k ) (35) 

where f s {k ) = co\{f s (j k )} € ft 6nj(fc) . Noting that 

f(k) = B(k)f s (k) (36) 

and using the principle of virtual work, it follows from Eq. (23) that the modal spatial forces 
f m (k ) G for the k th body are 



/.(*) = 


IP(*)/.(*) 

m 


Premultiplying Eq. (35) by ^ B(k) ^ an< ^ Eq. (25), Eq. (29), and Eq. (37) leads to the 

following recursive relationship for the modal spatial forces: 

f m (k) = $(k,k — 1 )fm{k — 1 ) + M m (k)a m (k) + b m (k) + K m (k) / d(k) ( 38 ) 


where 


bm(k) = 


A ( Yl*{k) 


[b(k) + M,(fc)a(fc)] G (k) 


and the modal stiffness matrix 


K r k \ £ ( n*(fc)A' s (fc)n(fc) o\ € ^ (fc )xA7(*) 


The expression for K m (k) in Eq. (40) uses the fact that the columns of B*(k) are the deformation 
dependent rigid body modes for the k th body, and hence they do not contribute to its elastic strain 
energy, ^^hen a deformation dependent structural stiffness matrix /\ 5 (/j) is used, 

K s (k)B*(k) = 0 (41) 

However, common practice (followed in this paper) uses a constant, deformation-independent struc- 
tural stiffness matrix. This leads to the apparently anomalous situation wherein Eq. (41) does not 
hold exactly. All these fictitious extra terms on the left-hand side of Eq. (41) are commonly ignored. 
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The velocity-dependent bias term b m (k) is formed using modal integrals generated by standard 
finite-element programs, and a detailed expression for it is given in [6]. From Eq. (38), the operator 

expression for the modal spatial forces f m = col|/ m (fc)| E ft^ for all the bodies in the chain is 

fm = HM m a m + b m + K m ti) (42) 

where 

M m = diag{M m (fc)} € K m = diag{tf m (*)} G , and b m = col {b m (k)} € ft* 

From the principle of virtual work, the generalized forces vector T G ft* for the multibody system 
is 

T = Ufm (43) 

2.5 Operator Expression for the System Mass Matrix 

Collection of the operator expressions in Eq. (22), Eq. (28), Eq. (42) and Eq. (43) leads to: 

Vm = 

a m = **(W** + ««) (44) 

fm = $(M m a m + b m + K m d) = + $(M m f *a m + b m + K m d) 

T - Hfm = n$M m $-H*x + H$(M m **a m + bm) 

= Mx + C 

where 

M = H$Mm**H* G and C = ?f$(M m $*a m + b m + K m d) € ft* (45) 

Here M is the system mass matrix. The expression 7i* is referred to as the Newton-Euler 

Operator Factorization of the mass matrix. The term C is the vector of Coriolis, centrifugal, and 
elastic forces for the system. 

The operator expressions for M and C are identical in form to those for rigid multibody systems 
(see [1,7]). This similarity is extremely useful in the extension of recursive algorithms from rigid 
multibody systems to flexible multibody systems. 


3. Composite Body Forward Dynamics Algorithm 

The forward dynamics problem for a multibody system requires computing the generalized accel- 
erations x for a given vector of generalized forces T and state of the system {i?,x}- The composite 
body forward dynamics algorithm described below consists of (a) computing the system mass matrix 
M, (b) computing the bias vector C, and (c) solving the linear matrix equation for x: 

Mx = T-C (46) 

Section 4 describes the recursive articulated body forward dynamics algorithm that does not require 
the explicit computation of either M or C. 


Lemma 3.1: Define the composite body inertias R(k) G ^ A ' r ( i ) xAA ( fc ) recursively for all the 

bodies in the serial chain as follows: 


r R( 0) = 0 

for k = 1 • • • N 

' R(k ) = $(k,k-l)R(k-l)$‘(k,k-l) + M m (k) 

end loop 


(47) 
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Also define R = diag{i£(fc)} € ^ x77 - Then, observe the following spatial operator decomposition 
where $ = $ — /: 

= R + i>R + R$* ( 48 ) 

Physically, R(k) is the modal mass matrix of the composite body formed from all the bodies 
outboard of the k th hinge by freezing all their (deformation plus hinge) degrees of freedom. It 
follows from Eq. (45) and Lemma 3.1 that 

M = HSM m W = HRH' + nSRH* + HRi>*H* (49) 

Note that the three terms on the right of Eq. (49) are block diagonal, block lower triangular and 
block upper triangular respectively. The algorithm for computing the mass matrix M computes 
these terms recursively. The main recursion proceeds from tip to base, and computes the blocks 
along the diagonal of M. As each such diagonal element is computed, a new recursion to com- 
pute the off-diagonal elements is spawned. Its structure is similar to that of the composite body 
algorithm for computing the mass matrix of rigid multibody systems (see [8,9]), and is as follows: 

[ R{ 0) = 0 

I for k = 1 ■ • • N 

$(/:,& — 1 )R(k — 1 )$*(&,& — 1) + M m (k) 

A(k)<t>(t k ,k - l)R TT {k- l)<f> m (t k ,k - l)A*(k) + M m (k) 

R(i k)H*(k) 

H{k)X(k) 

r 

for j = (k + 1) • • - N 

< X(j) = *(j,j-l)X(j-l) = A(mt j J-l)X r (j-l) 

M(j, k) = M m (k, j) = H(j)X{j) 

end loop 

(50) 

The structure of the above algorithm for computing the mass matrix closely resembles the composite 
rigid body algorithm for computing the mass matrix of rigid multibody systems [8,9]. Like the 
latter, it is also highly efficient. Additional computational simplifications of the algorithm arising 
from the sparsity of both W/(fc) and H r {k) are easy to incorporate. 


m = 

X(k) = 
< M s {k,k) = 

end loop 


4. Factorization and Inversion of the Mass Matrix 

An operator factorization of the system mass matrix M, referred to as the Innovations Operator 
Factorization , is derived. This factorization is an alternative to the Newton-Euler factorization in 
Eq. (45). In contrast with the latter, the factors in the Innovations factorization are square and 
invertible. Operator expressions for the inverse of these factors lead to an operator expression for 
the inverse of the mass matrix. Use of further operator identities results in the recursive articulated 
body forward dynamics algorithm in Section 5. The operator factorization and inversion results 
here closely resemble those for rigid multibody systems (see [1]). 

The following recursive algorithm defines some required articulated body quantities. This algorithm 
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has the structure of the Riccati equation of Kalman filtering theory [9]: 


' P+( 0) = 

for k = 1 •••N 
P(k) = 
D(k) = 
G(k) = 
< K(k+l,k) = 
T(k) = 
P+(k) = 
V(k+l,k) = 
. end loop 


0 


$(k,k- l)P+(k- l) + Mm(k) G ^)xW(k) 

H(k)P{k)H*{k) G 8^< fc > x ^W 
P{k)H*(k)D- l {k ) G 
$(*-H, *)<?(*) 

7 - G(k)H(k ) G 
r(fc)P(ifc) 6 &(*)>&{*) 

$(k + i,k)T(k)e$Fw xjr w 


(51) 


The operator P G JJ^ X ^ is defined as the block diagonal matrix with the k tfl diagonal element 
being P(k). The quantities defined in Eq. (51) form the component elements of the following spatial 
operators: 


D 

A 

HPH* = diag{7)(fc)} G ^ xN 

G 

A 

PH'D~ l = diag{(7(fc)} G 

Ii 

A 

4Ge $F xAr 

T 

A 

I - GH - diag{r(*)} G %F xJf 

£<$ 

A 

£ 9 t G $P xT 


(52) 


The only nonzero block elements of K and £$ are the elements K (k + 1, k) and + 1, k) respec- 
tively along the first sub-diagonal. The spatial operator G is formed by a set of spatial Kalman 
gains [9] for a spatially recursive Kalman filter. 

As in the case for is nilpotent, so the operator $ can be defined as 


( 

* = (7-&J,) -1 = 


7 

*( 2 , 1 ) 


0 

7 


n 


e $F xjr 


V¥(JV, 1) ¥(A,2) ... I) 


(53) 


where 

^(i,j) =¥(*,*- 1) ••• ¥(j + l,j) for i > j 

The structure of the operators Sq and $ is identical to that of the operators £$ and § respectively 
except that the component elements are now ^(ij) rather than 4>(z, j). Also, the elements of 
$ have the same semigroup properties as the elements of the operator $, and as a consequence, 
high-level operator expressions involving them can be directly mapped into recursive algorithms, 
and the explicit computation of the elements of the operator $ is not required. 

The Innovations Operator Factorization of the mass matrix is defined in the following lemmas 
established in [1]. 


Lemma 4.1: 

M = [I + H$K]D[I + (54) 
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Note that the factor [/ + H*K] 6 is square, block lower triangular and nonsingular, while 

D is a block diagonal matrix. This factorization may be regarded as a block LDL’ decomposition 
of M. The following lemma gives the closed form operator expression for the inverse of the lactor 

[/ + H*K). 


Lemma 4.2: 

[/ + HQhT' = [/ - H9K] ( 55 ) 


Lemma 4.3: 

M-'=[I-H9KYD- l [I-H9K) (56) 

Once again, note that the factor [/ - H9K] is square, block lower triangular and nonsingular 
and so Lemma 4.3 may be regarded as providing a block LDL * decomposition of M . This 
decomposition however is model-based, in the sense that the physical model of the system is used 
to conduct computations. This means that every step in the decomposition has a corresponding 
physical interpretation which adds a substantial amount of insight into the decomposition. 


5. Articulated Body Forward Dynamics Algorithm 

The operator- based mass matrix inverse leads to a recursive forward dynamics algorithm. The 
structure of this algorithm is completely identical in form to the articulated body algorithm for serial 
rigid multibody systems. Its structure is that of a Kalman filter and a Bryson-Fraz.er smoother [9j. 

The following lemma, established in [1], describes the operator expression for the generalized ac- 
celerations x in terms of the generalized forces T . 


Lemma 5.1: 

x = [l - Ti^K) m D~ x [t - HV{KT + Pa m + b m + A' m t?}] - A (57) 

As in the case of rigid multibody systems [1,10], the direct recursive implementation of Eq. (57) 
leads to the following recursive forward dynamics algorithm: 

z + (0) = 0 

for k = 1 • • • n 

*(*) = *(fc,fc - l)z + {k - 1) + P(k)a m (k) + b m (k) + Km(k)d(k) 

< c(k) = T(k)-n(k)z(k) 
y{k) = D~'(k)i(k) 
z + (k ) = z(k) + G(k)<(k) 

end loop 

(58) 

a m (n + 1 ) = 0 

for k = n • • 1 

a +(fc) - * m (k+Uk)a m (k+l) 

* x(k) = v(k) - G*(fc)a+(fc) 

a rn ( k ) = ttj,(l.') + 'H“(k)x(k) + a m (k) 

end loop 

All the degrees of freedom for each body are characterized by its joint map matrix H'{.) and ar p 
processed together at each recursion step in this algorithm. However, by taking advantage of the 
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sparsity and special structure of the joint map matrix, additional reduction in computational cost 
is obtained by processing the flexible dofs and the hinge degrees of freedom separately. These 
simplifications are described in the following sections. 

Instead of giving detail, the conceptual approach to separating modal and hinge degrees of freedom 
is described. First, recall the velocity recursion equation in Eq. (13) 

V m (k) = *'(k + 1 % k)V m (k + 1) + H*(k) X (k) (59) 

and the partitioned form of H(k) in Eq. (15) 


H(k) = 


«/(*) ^ 

Hr{k) ) 


( 60 ) 


Introducing a dummy variable rewrite Eq. (59) as 

V m (k') = ^(ir + + 1) + 

v m (k) = **(k',k)v m (k') + n;(k) 0 (k) (6i) 


where 

${k + 1 ,k ( ) = $(k + 1 ,k) and <!>(&', fc) = I 

Conceptually, each flexible body is now associated with two bodies. The first one has the same 
kinematical and mass/inertia properties as the real body and is associated with the flexible degrees 
of freedom. The second body is a fictitious body, and it is massless and has zero extent. It is 
associated with the hinge degrees of freedom. The serial chain now contains twice the number of 
bodies as the original one, with half the new bodies being fictitious. The new H* operator has 
the same number of columns but twice the number of rows as the original H* operator. The new 
$ operator has twice the number of rows as well as twice the number of columns as the original. 
An analysis similar to those of the previous sections leads to an operator expression similar to 
Eq. (57). This implies a recursive forward dynamics algorithm like Eq. (58). However each sweep 
in the algorithm now contains twice as many steps as the original algorithm. But since each step 
now processes a smaller number of degrees of freedom, the overall computational cost is reduced. 


5.1 Simplified Articulated Body Forward Dynamics Algorithm 

The complete recursive articulated body forward dynamics algorithm for a serial flexible multibody 
system follows from recursive implementation of Eq. (57). The algorithm has the following steps: 
(a) compute the articulated body quantities, (b) do a base-to-tip recursion for the modal spatial 
velocities V m (k) and the bias terms a m (k), b m (k ), and (c) do a tip-to-base recursion followed by a 
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base-to-tip recursion for the joint accelerations 


' 

p£( 0) = o 

for k = 

1 •■■N 

r(*) 

1 

* 

1 

1 

II 

p(k) 

= A(k)T(k)A'(k) + M m (k) 

D,(k) 

= Hj(k)P(k)Hj(k) 

n(k) 

= [P r Hk), P rr (k)]H’j(k) 

9(k) 

= n(k)Dj\k) 

PR(k) 

= P”(k)-g(k)n'(k) 

D R (k) 

= Hr{k)P R {k)Hy{k) 

G R (k) 

= P R {k)Hy(k)D- R \k ) 

rn(k) 

= I - G R (k)IIjr(k) 

p£(k) 

end loop 

= r R (k)P R (k) 


( 62 ) 


for k 


4(0) = o 

1 •■•N 


*(k) = 


cj(k) = 

vj{k) = 

*n{k) = 

<n(k) = 

V R {k) = 

4 (*) = 

end loop 


f */(*> 1 

( 'r(k) ) 

A(k)<t>(h.,k - 1)4(* - 1) + M*) + A ’ m (k)d(k) 6 
Tj(k ) - Zj(k) + [n^(fc)]‘*r(fc) € % n '" [k] 
j l (k)( / (k)ex n * 


D~ i 


-(*) 


*(*) + $(**/(*) + PR(k)a mR (k) e 5? 

Tfi(l-) - Hr(k)z R (k) e ^ ( *> 

*n(*) + € » 6 


a, 


for k 

G^( ) 

/*(*) 

*K*) 

m ( ) 

end loop 


N' 


,(* + 1) = 0 

= <t> m (tk+ \,k)A*{k + l)a m (k + 1) € 
= »R(k)-G* R (k)a+(k)6ft n 'W 
= a R (k ) + H)r(kmk) + a mfl (*) G 


(63) 


The recursion in Eq. (63) is obtained by carrying out simplifications of the recursions in Eq. (58) 
in the same manner as described in the previous section for the articulated body quantities. 

In contrast with the composite body forward dynamics algorithm of Section 3, this algorithm 
docs not explicitly compute either M or C. This algorithm is similar to those for rigid multibody 
systems [1,11]. 
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6. Computational Cost 

The computational costs of the composite body and the articulated body forward dynamics algo- 
rithms are compared. For low-spin multibody systems, it has been suggested in [12] that using 
ruthlessly linearized models for each flexible body can lead to significant computational reduction 
without sacrificing fidelity. These linearized models are considerably less complex than the lull 
nonlinear models and do not require much of the data on modal integrals for the individual flexible 
bodies. All computational costs given below are based on the use of ruthlessly linearized models 
and the computationally simplified steps described in [6]. 

6.1 Computational Cost of the Composite Body Forward Dynamics Algorithm 

The composite body forward dynamics algorithm described in Section 3 is based on solving the 
linear matrix equation 

M\ = T-C 

The computational cost of this forward dynamics algorithm is as follows: 


1. The cost of computing R{k) for the k th body by using the algorithm in Eq. (50) is 
[48n m (/c) + 90] M + [n? u (fc) + f n m (k) + 116]A. 


2. The contribution of the k ,h body to the cost of computing M (excluding cost of R(k) s) using 
the algorithm in Eq. (50) is 

{k\V2n 2 m (k) + 34n m (fc) + 13]} M + {k[[ln 2 m (k) + 24n m (fc) + 13]} A. 


3 Setting the generalized accelerations \ = 0, the vector C can be obtained by using an inverse 
dynamics algorithm for computing the generalized forces T. The contribution of the k body 
to the computational cost for£(A:) is {2 n^(k) + 54n m (fc) + 200} M + {2n m (fc) + 50 n m (k) + • } 


4. 

5. 


The cost of computing T - C is {Af} A. 


The cost of solving the linear equation in Eq. (46) for the accelerations \ is 

{IA' 3 + \\ ri - |A'} M + {iA' 3 + A' 2 - | A'} A. 


The overall complexity of the composite body forward dynamics algorithm is 0( A" 3 ). 


6.2 Computational Cost of the Articulated Body Forward Dynamics Algorithm 

The articulated body forward dynamics algorithm is based on the recursions described in Eq. (62) 
and Eq. (63). Since the computations in Eq. (47) can be done prior to the dynamics simulation, 
the cost of this recursion is not included in the cost of the overall forward dynamics algorithm 
described below: 


1 The algorithm for the computation of the articulated body quantities is given in Eq. (62). The 
' step involving the computation of D~ l (k) can be carried out either by an explicit inversion of 
D(k) with Oin^k)) cost, or by the indirect procedure described in Eq. (62) with U(n m (k)) 
cost. The first method is more efficient than the second one for n m {k) < 7. 


• Cost of Eq. (62) for the k th body, based on the explicit inversion of D(k) (used when 

n m {k) < 7), is . 

{§ n 3 m (k) + f n 2 m (k) + &n m (k) + 18o} M + {fn^tfc) + + ^n m (fc) + 161) A. 

• Cost of Eq. (62) for the k th body based on the indirect computation of D~ l (k) (used 
when n m (k) > 8) is { 12 n^(fc) + 255n m (fc) + 572} M + { 13 n 2 n (k) + 182 n m (k) + 4 15} A. 
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2. The cost for the tip-to-base recursion sweep in Eq. (63) for the k th body is 
{n 2 m {k) + 25 n m (k) + 49} M + {n^(fc) + 24 n m (k) + 50} A. 

3. The cost for the base-to-tip recursion sweep in Eq. (63) for the k ,h body is {18 n m (k) + 52} M + 
{I9n m (k) + 42} A. 

The overall complexity of this algorithm is 0(An^,), where n m is an upper bound on the number 
of modes per body in the system. 

The articulated body algorithm is more efficient than the composite body algorithm as the number 
of modes and bodies in the multibody system increases. Figure 2 contains a plot of the computa- 
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Figure 2: A comparison of computational costs for the forward dy- 
namics algorithms for a flexible multibody serial chain 
system with 10 flexible bodies. 

tional cost (in floating point operations) versus the number of modes per body for a serial chain 
with ten flexible bodies. The articulated body algorithm is faster by over a factor of 3 for 5 modes 
per body, and by over a factor of 7 for 10 modes per body. The divergence between the costs for 
the two algorithms becomes even more rapid as the number of bodies is increased. 


7. Extensions to General Topology Flexible Multibody Systems 

Extension to general tree and closed-chain systems is similar to methods given in prior results for 
rigid body configurations [7]. The key is that the operator description does not change as the 
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topology changes. Extending the serial chain results of this paper to tree topologies takes the 
following steps: 

1. For any outward base to tip(s) recursion, at each body, the outward recursion must be con- 
tinued along each outgoing branch emanating from the current body. 

2. For an inward tip(s) to base recursion, at each body, the recursion must be continued inward 
only after summing up contributions from each of the incoming branches of the body. 


A closed-chain flexible multibody system can be regarded as a tree topology system with additional 
closure constraints [7]. 

8. Conclusions 

This paper uses spatial operator methods to develop a new dynamics formulation and spatially 
recursive algorithms for flexible multibody systems. The operator description of the flexible system 
dynamics is identical in form to the corresponding operator description of the dynamics of rigid 
multibody systems. A significant advantage of this unified approach is that it allows ideas and 
techniques for rigid multibody systems to be easily applied to flexible multibody systems. All 
of the computations are mechanized within a spatially recursive Kalman filtering and smoothing 
architecture. An extension of this algorithm to handle prescribed motion is described in reference 
[13]. 

The computational efficiency of the dynamics algorithms described in this paper makes it possible to 
implement real-time, high-fidelity, hardware-in-the-loop simulation of complex multibody systems 
such as spacecraft, robot manipulators, vehicles etc. Such simulations are essential during the 
design and testing of control and fault recovery algorithms. The articulated body forward dynamics 
algorithm is currently being used to simulate the dynamics of planetary spacecraft. One application 
is a spacecraft currently being assembled for a comet and asteroid rendezvous mission [14]. The 
multibody model for the spacecraft is of tree topology, and consists of a flexible central bus with 9 
articulated appendages and 22 hinges’ degrees of freedom. The simulation software provides a new 
capability for high speed simulation of the spacecraft. A real-time version has also been developed. 
Validation of this software was carried out by running independent simulations of the spacecraft 
using a standard flexible multibody simulation package [15]. Results from the two independent 
simulations show complete agreement. 
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